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INTRODUCTION 


In  a  previous  report, ^  this  writer  presented  a  finite  element-variational 
formulation  which  discretizes  the  spatial  and  time  variable  in  the  same 
manner.  The  method  was  applied  to  a  problem  of  beam  motion  subjected  to 
moving  concentrated  forces.  Results  were  shown  to  be  in  excellent;  agreement 
with  known  solutions.  This  same  formulation  is  now  applied  to  the  problem  of 
a  moving  couple,  i.e.,  a  concentrated  bending  moment. 

A  recent  investigation  by  S.  H.  Chu^  on  the  interacting  forces  between  a 
projectile  and  the  cannon  tube  indicates  that  the  couple  produced  by  the 
eccentricity  of  the  projectile  as  it  moves  down  the  tube  may  be  of  such  a 
magnitude  that  its  effect  on  the  tube  motion  becomes  significant.  It  is  then 
important  that  the  problem  associated  with  moving  moments  can  be  analyzed 
adequately.  The  purpose  of  this  note  is  to  present  the  modification  necessary 
to  the  previous  formulations  so  that  the  solutions  of  a  beam  motion  problem 
under  a  moving  bending  moment  can  be  obtained  routinely.  Results  of  a 
cantilevered  beam  subjected  to  such  a  load  are  also  present. 

DIFFERENTIAL  EQUATION  AND  NONDIMENSIONALIZATION 

Consider  a  Euler-Bernoulli  beam  subjected  to  a  moving  couple  M.  The 
equation  differential  can  be  written  as 

^J.  J.  Wu,  ’‘Beam  Motions  Under  Moving  Loads  Solved  by  Finite  Element  Method 
Consistent  in  Spatial  and  Time  Coordinates,”  ARLCB-TR-80046 ,  USA  ARRADCOM, 
Large  Caliber  Weapon  Systems  Laboratory,  Benet  Weapons  Laboratory, 

Watervliet,  NY,  November  1980. 

^S.  H.  Chu,  ”In-Bore  Motion  Analysis  of  the  155  mm  XM712  Projectile  When  Fired 
in  the  M198  Howitzer,”  Proceedings  of  the  Army  Symposium  on  Solid  Mechanics, 
AMMRC  MS  80-4,  Army  Materials  and  Mechanics  Research  Center,  Watertown,  MA, 
pp.  270-288,  1980. 
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Ely***'  +  pAy  =  -M6’(x-x)  .*  (1) 

where  y(x,t)  denotes  the  beam  deflection  as  a  function  of  spatial  coordinate  x 
and  time  t,  E,  I,  A,  p  denote  elastic  modulus,  second  moment  of  inertia  area 
and  material  density  respectively.  A  dirac  function  is  denoted  by  6,  x  =  x(t) 
is  the  location  of  M,  a  prime  (’)  denotes  differentiation  with  respect  to  x 
and  a  dot  (*),  differentiation  with  respect  to  t.  Note  that  the  right  hand 
side  of  Eq.  (1)  has  a  dimension  of  force  due  to  the  fact  that 

d  -  - 

6’(x-x)  =  —  [6(x-x)] 
dx 

and  it  has  a  dimension  of  (length)“^. 

Introducing  nondimensional  quantities 

^  ^  ^ 

y  «  y/i  ,  X  =  x/Ji  ,  t  (2) 

where  I  is  the  length  of  the  beam  and  T  is  a  finite  time,  within  0  <  t  <  T, 
the  problem  is  of  interest,  Eq.  (1)  can  be  written  as 


^  y2y  -  _  q6’(x-x) 

The  hats  ('")  have  been  omitted  in  Eq.  (3)  and 

c 

y  =  “ 

T 


with 


Q  = 

El 

pA£^ 


El 


(3) 


(4) 
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Boundary  conditions  associated  with  Eqs,  (1)  or  (2)  will  now  be  introduced  in 
conjunction  of  a  variational  problem.  Consider 

61  =  0  (5a) 


with 


11  •  •  —  — 

I  =  /^[y’'y*"  -  Y^yy*  +  Q6(x-x)y*]dxdt 


+  dt{kiy(0 ,t)y*(0 ,t)  +  k2y’ (0 , t) y*' (0 ,t) 

+  k3(y(l ,t)y*(l ,t)  +  k4y' ( 1 ,t)y* ' ( 1 ,t) } 

+  dx{k5[y(x,0)  -  Y(x) ] y*(x,l) }  (5b) 

where  y*(x,t)  is  the  adjoint  variable  of  y(x,t).  If  one  takes  the  first 
variation  of  I  considering  y(x,t)  to  be  fixed: 

(6l)5y,o  “0  (5a') 

and  consider  6y*  to  be  completely  arbitrary,  it  is  easy  to  see  that  Eqs.  (5) 
is  equivalent  to  the  differential  Eq.  (3)  and  the  following  boundary  and 
initial  conditions. 

y"'(0,t)  +  kiy(0,t)  =  0 
y"(0,t)  -  k2y'(0,t)  =  0 

0  <  t  <  1  (6a) 

y"'(l,t)  -  k3y(l,t)  =  0 
y"(l,t)  +  k4y'(l,t)  =  0 

. 

y(x.O)  =  0 

and  .  0  <  X  <  1  (6b) 

y(x,l)  -  k5[y(x,0)  -  Y(x)]  =  0 
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Taking  appropriate  values  for  ,  k2,  k3,  and  k4,  problems  with  a  wide  range 
of  boundary  conditions  can  be  realized.  The  initial  conditions  in  Eqs.  (6b) 
are  that  the  beam  has  zero  initial  velocity,  and,  if  one  takes  k5  to  be  »  (or 
larger  number  compared  with  unity) , 

y(x,0)  «  Y(x) 

The  meaning  for  cases  where  k5  is  not  so  need  not  be  our  concern  here* 

To  derive  the  finite  element  matrix  equations,  one  begins  with  Eq.  (5a’) 


and  write 


('SDSy-O  *  0 


(7a) 


+  k3y(l ,t)6y*(l ,t)  +  k4y' ( 1 , t) 5y* ' ( 1 , t) ] 


(7b) 


Introducing  element  local  variables 


(i) 

5  =  C  =  Kx-i+1 


(i) 

n  =  n  =  Lt-j+i 


(8a) 


or 


1 

X  =  -  (5+i-l) 


K 


(8b) 


1 


t  =  -  (5+j-l) 


L 
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where  K  is  the  number  of  divisions  in  x  and  L,  in  t.  (A  typical  grid  scheme 
is  shown  in  Figure  1).  Equation  (7b)  can  now  be  written  as 
K  L  1  1  y^L 

I  I  /  /  [~  y(ij)«y*(ij)]‘i5dTi 

isal  U  U  L  K 


^  ,1  M 

+  I  /  dn  [—  y(ij)(0,n)6y*(ij)(0,Ti)  +  k2  —  y'(ij)(0,n)6y*'(ij)(0,n) 
j=sl  0  L.  L 

+  I  /g  --  [ir^k5(y(ij)(C,0)6y*(ij)(5,l))] 


K  L  Q  1  1  _ 

=  -  I  I  -  L  L  6'(x-x)6y*(ij)(5,n)d5dn 

i*l  j-l  L  0  0 


K  Y^ks  1 

+  1  -  /„  dC  [Y(i)(5)6y*(iL)(5,l)]  ,  (9) 

i»l  K  0 

The  shape  function  vector  is  now  introduced.  Let 

y(ij)(C,n)  =  aT(5,n)Y(ij) 


y*(ij)(5,Tl)  =  aT(5jn)Y*(ij)  =  Y*T(ij)a(5,n)  (10) 
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Equation  (9)  then  becomes 


K  L 


.K 


I  I  t--  A  -  ?)  !(ij) 


i=l  j=,l  -  ■  L,  ~  K 

L  ki  k2K^ 

+  I  {--  Bi  +  -- —  B2}  Y(ij) 

i=l  ~  L  ~  L  -  - 

L  k3  k^K^ 

+  I  5Y*'^(Kj)  B3  +  --  B4}  Y(ij) 

i=l  -  L  -  L  -  ~ 

K  T^k^ 

+  I  5Y*T(iL)  {---  B5}  Y(iL) 


i=l 


IC  Ij  ^ 

I  I  7  F(ij)  +  I  ^*'^(±L)  ““  5(i) 

i-l  j=l  -  L  1=1”  K- 


(11) 


where,  as  it  can  be  seen  easily,  that 


A  “  /  /  ^  E.E  EE 

0  0  ~ 


and 


.1  .1  _ 

B=  /  /  a^alv,  d^dn 

Bi  =  f  a(0  ,ri)a'B(0  ,ri)dn  ,  ^2  =  a^  ^(0 ,  n)a^^  ^(0  ,  n)dn 


B3  =  /  a(l  ,n)a'^(l  ,n)dn  ,  84=/  a  r(l  ,n)a^  r(l  ,n)dn 

-n~  ~  -  Dr’ 


,1 


(12) 


B5  =  a(5,l)aT(5,0)d5 


1  1  -  -  ,1 

F(ij)  =  /q  /q  a(5,n)(S'(ij)(5-C)d5dTi  ,  G(i)  -  a(5,l)Y(i)(C)d5 
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where 


d  - 


«'(lj)(M)  -  --  «(ij)(?-5) 


is  the  local  version  of  the  function  6(x-x)  appeared  in  Eq.  (9).  The  specific 
form  of  6(y)(5-5)  will  be  described  later  in  a  paragraph  prior  to  Eqs,  (18). 
Now  Eq.  (11)  can  be  assembled  in  a  global  matrix  equation 


6Y*T  K  Y  =  F 


(13) 


By  virtue  of  the  fact  that  6Y*  is  not  subjected  to  any  constrained  conditions, 
one  has  ^ 


(14) 


K  Y  =  F 


which  can  be  solved  routinely.  Numerical  results^  of  several  problems  in  this 
class  will  be  presented  in  a  later  section. 

< 

FORCE  VECTOR  DUE  TO  A  MOVING  COUPLE 

We  shall  describe  here  the  procedures  involved  to  arrive  at  the  force 
vector  contributed  by  a  moving  couple.  This  force  vector  has  appeared  in  Eq. 
(12)  as 


(15a) 


Perform  integration-by-parts  once.  Equation  (15a)  can  be  written  as 


F(ij)  =  -  /q  ^0 


(15b) 
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The  shape  function  a(5>n)  is  a  vector  of  16  in  dimension.  In  the  present  for¬ 
mulation  we  have  chosen  the  form: 

=  bi(5)bj(n)  ,  k  =  1,2,3, ...16  (16a) 

i,j  =«  1,2, 3, 4 

and 

=  bi’(5)bj(5)  (16b) 

The  relations  between  k  and  i,j  are  given  in  Table  I.  These  are  the  conse¬ 
quences  of  the  choice  of  the  shapq  function  $uch  that  the  generalized 

coordinates  of  the  (ij)th  element,  represent  the  displacement,  slope,  veloc¬ 
ity,  and  angular  velocity  at  the  local  nodal  points.  Thus 

ti(C)  =  I  ;  bi’(5)  =  I  (17) 

p=l  p=l 

The  values  of  bj^p  are  given  in  Tables  II  and  III, 


TABLE  I.  EELATIONSHIP  BETWEEN  (i,j)  AND  k  IN  EQUATION  (16) 


k 

(i,j) 

k 

(i,j) 

1 

(1,1) 

9 

(1,3) 

2 

(2,1) 

10 

(2,3) 

3 

(1,2) 

11 

(1,4) 

4 

(2,2) 

12 

(2,4) 

5 

(3,1) 

13 

(3,3) 

6 

(4,1) 

14 

(4,3) 

7 

(3,2) 

15 

(3,4) 

8 

(4,2) 

16 

(4,4) 
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TABLE  II.  VALUES  OF  bip  IN  EQUATION  (17) 


p 

i 

1 

2 

3 

4 

1 

1 

0 

-3 

1 

2 

0 

1 

-2 

1 

3 

0 

0 

3 

-2 

4 

0 

0 

-1 

1 

TABLE  III.  VALUES  OF  b'ip  IN  EQUATION  (17) 


Now,  let  us  consider  6^^j^(5“5)*  This  function  represents  the  effect 
of  the  Dirac  delta  function  6(x-x)  on  the  (ij)th  element.  If  the  curve  of 
travel  x  =  x(t)  does  not  go  through  the  element  (i,j) ,  ^(ij)(5“5)  ~  0.  If  If 
passes  through  that  element ,  one  has 

6(ij)(5-i)  =  S(x-x)  =  K6(5-5)  (18a) 
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with 


5  =  5(n) 

The  function  C(n)  is  derived  from  x  =  x(t).  For  example,  if  the  force 
with  a  constant  velocity,  one  has 


X  =  x(t)  =  vt 


it  follows  from  Eqs.  (8)  that 


vK 


5  =  5(n)  =*  -i+1  +  —  (h+j-1) 

L 

With  Eqs.  (16),  (17),  (18),  and  (19),  one  writes  (15)  as 


,1  ,l 


^(ij)k  ^  a]^,5(5  ,n)'5(5-C)d5dn 


,1  ,1  -  -  p-1  q-1-  - 

F(ij)k  =  K/q  /q  b'pbjq  C  n  6(?-5)d5dn 

Equation  (20)  can  then  be  evaluated  easily  once  the  exact  form  of  C  is 
written.  For  example,  if  5  =  h,  Eq.  (20)  reduces  to 


4  4 


^(ij)k  “  I  I  ^’ip  ^iq  /  5 

p=l  q=l 


1  _p*-q-2 


dK 


^44  kb’ipbjq 

"  pil  qii  "^q=r" 


(18b) 

moves 

(19a) 

(19b) 

(20a) 

(20b) 

(21) 
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TABLE  IV.  DEFLECTION  y(x,t)/A  OF  A  CANTILEVERED  BEAM  UNDER  A 
MOVING  CONCENTRATED  MOMENT  (T  =  10^°  sec.) 


x/i 

t/T 

0 

0.25 

0.50 

0.75 

1.00 

0. 

0. 

0. 

0. 

0. 

0. 

0.25 

0. 

.03125 

.09375 

0.15625 

0.21875 

0.50 

0. 

.03125 

.12500 

0.25000 

0.37500 

0.75 

0. 

.03125 

.12500 

0.28125 

0.46875 

1.00 

0. 

.03125 

.12500 

0.28125 

0.50000 

TABLE  V.  DEFLECTION  y'(x,t)/;i  OF  A  CANTILEVERED  BEAM  UNDER  A 
MOVING  CONCENTRATED  MOMENT  (T  =  10^°  sec.) 


x/Z 

t/T 

0 

0.25 

0.50 

0.75 

1.00 

0. 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.25 

(0.0072) 

0.2366 

0.2481 

0.2496 

0.2500 

0.50 

(0.0021) 

0.2481 

0.4856 

0.4981 

0.5000 

0.75 

0.0005 

0.2496 

0.4981 

0.7366 

0.7500 

1.00 

0.0000 

0.2505 

0.5021 

0.7572 

1.0000 
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TABLE  VI.  DEFLECTION  y(x,t)/i!.  OF  A  CANTILEVERED  BEAM  UNDER  A  MOVING  CONCENTRATED  MOMENT 


o 

0) 

Cfi 
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TABLE  VII.  DEFLECTION  y(x,t)/Jl  OF  A  CANTILEVERED  BEAM  UNDER  A  MOVING  CONCENTRATED  MOMENT 
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TABLE  VIII.  DEFLECTION  y(x,t)/Jl  OF  A  CANTILEVERED  BEAM  UNDER  A  MOVING  CONCENTRATED  MOMENT 
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NUMERICAL  DEMONSTRATIONS 

Some  numerical  results  obtained  will  now  be  presented.  Let  us  consider  a 

cantilevered  beam  subjected  to  a  unit  moving  couple  with  a  constant  velocity 

£ 

V  =  - 

T 


As  T  varies  from  "  to  0,  the  velocity  varies  from  0  to 

It  will  be  helpful  to  compare  v  with  some  reference  velocity  which  is  a 


characteristic  of  the  given  beam.  It  is  known  that  for  a  cantilevered  beam, 
the  first  mode  of  vibration  has  a  frequency  (see,  for  example,  reference  3) 
u)  1  ,1.875^,  0*560 

fi  =  —  =  —  ( - J  = -  (cycles  per  seconds) 

2ti  2ir  c  c 


and  the  period, 


Ti  =  1.786  c 


El 


Consider  the  vibration  as  standing  waves.  They  travel  at  a  speed 


1.12£ 

vi  =  2£fi  = - — 

c 


Hence,  the  relative  velocity 

-  V  c 

V  =  —  =  —  =  0.893  - 

v^  2T  T 


Or.  N.  Tong,  Theory  of  Mechanical  Vibration,  John  Wiley,  New  York,  1960,  p. 
257;  p.  256. 
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We  shall  take  c  »  1.0  for  the  moving  force  problems.  Thus,  »  0.560  Hz. 

=  1.786  sec.  and 

V  =*  0.893/T 

Using  grid  schemes  of  4x4  (i.e.,  four  segments  in  spatial  and  four  in 
time  coordinates)  and  8x4.  Tables  IV  through  VIII  show  the  beam  deflections 
(and  slopes)  as  the  concentrated  moment  Q  =  1.0  moves  from  the  left  to  the 
right  end.  Since  we  have  defined  T  as  the  time  required  for  fhe  load  to 
travel  from  one  end  to  another,  t  “  0.5X,  for  example,  indicates  the  point  \in 
time  when  the  load  is  at  the  midspan  of  the  beam. 

In  Tables  IV  and  V,  T  is  set  to  10^®  sec.  which  is  extremely  large 
compared  with  the  beam  characteristic  time  of  T^  =  1.786  sec.  The  solution 
should  reduce  to  the  static  problem.  This  is  certainly  the  case  as  shown  in 
these  two  tables.  These  results  are  obtained  using  a  grid  scheme  of  4x4. 

For  results  shown  in  Tables  VI  through  VIII  an  8x4  grid  scheme  has  been 
used.  The  beam  deflections  for  T  -  10,  1.0,  and  0.1  seconds  are  shown  in 
Table  VI,  VII,  and  VIII  respectively. 

Finally,  these  deflection  curves  are  also  plotted  in  Figures  2  through 
10.  From  these  figures  and  the  tabulated  results,  one  observes  that  while 
some  of  the  results  are  extremely  good,  others  are  changing  so  rapidly  with 
respect  to  time  or  space  variable  that  an  assessment  on  their  accuracy  is  very 
difficult.  Hence,  further  investigations  on  numerical  convergence  of  these 
data  is  necessary. 
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Figure  1 


j 


A  Typical  Finite  Element  Grid  Scheme  Showing  the  (i,j)^^ 
Element  and  the  Global>  Local  Coordinates. 
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Figure  2. 


Deflection(s)  of  a  Cantilevered  Beam 
Under  a  Moving  Couple. 
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Figure  3.  Deflectlon(s)  of  a  Cantilevered  Beam  Under  a 
Moving  Couple. 
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Figure  4«  Deflection(s)  of  a  Cantilevered  Beam  Under  a 
Moving  Couple. 
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Figure  5.  Deflection(s)  of  a  Cantilevered  Beam  Under  a 
Moving  Couple. 
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Figure  7. 


Deflection (s)  of  a  Cantilevered  Beam  Under 
a  Moving  Couple, 
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Figure  8. 


Deflection (s)  of  a  Cantilevered  Beam  Under  a 
Moving  Couple. 
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Figure  9.  Deflection(s)  of  a  Cantilevered  Beam  Under  a 
Moving  Couple. 
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Figure  10.  Deflectlon(s)  of  a  Cantilevered  Beam  Under  a 
Moving  Couple. 
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